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ABSTRACT 

We investigate, using the spherical Jeans equation, self-gravitating dynamical equilibria sat- 
isfying a relation pjcr^ °^ which holds for simulated dark-matter haloes over their whole 
resolved radial range. Considering first the case of velocity isotropy, we find that this problem 
has only one solution for which the density profile is not truncated or otherwise unrealistic. 
This solution occurs only for a critical value of a^^^ - 35/18 - 1.94, which is consistent 
with the empirical value of 1 .9 + 0.05. We extend our analysis in two ways: first we introduce 
a parameter e to allow for a more general relation p/crf oc r^"; second we consider veloc- 
ity anisotropy parameterised by Binney's p{r) = 1 - cTg/cr^. If we assume /3 to be linearly 
related to the logarithmic density slope y(r) = -d In p/d In r, which is in agreement with sim- 
ulations, the problem remains analytically tractable and is equivalent to the simpler isotropic 
case: there exists only one physical solution, which occurs at a critical a value. Remarkably, 
this value of a and the density and velocity-dispersion profiles depend only on e and the value 
/3o = /3{r = 0), but not on the value /3oc = /3{r — > oo) (or, equivalently, the slope d/S/dj of 
the adopted linear /3-y relation). For e = 3, aa-n = 35/18 - 2/?o/9 and the resulting den- 
sity profile is fully analytic (as are the velocity dispersion and circular speed) with an inner 
cusp p oc f-O+iopo)/') a very smooth transition to a steeper outer power-law asymptote. 
These models are in excellent agreement with the density, velocity-dispersion and anisotropy 
profiles of simulated dark-matter haloes over their full resolved radial range. If e = 3 is a uni- 
versal constant, some scatter inySo ~ may account for some diversity in the density profiles, 
provided a relation p/crj oz r""'"' always holds. 
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1 INTRODUCTION 

It has long been recognised that A'-body studies of large-scale 
structure formation in cold dark matter (CDM) cosmologies pro- 
duce dark-matter haloes whose density profiles are remarkably 
similar in shape over a wide range of halo virial mass (e.g . 
Dubinski & Carlbers 1991; Crone et al. 1994; Navarro et al.'l99a 
1997; Moore et al. 1999; Bullock et al. 2001). This 'universal' halo 
density distribution is characterised by a relatively shallow power- 
law behaviour in the inner parts, p ~ r^'', with y x \ typically 
inferred at the smallest resolved radii, which steepens gradually to 
an extrapolated y a; 3 - 4 at arbitrarily large radii. 

A physical explanation, based on first principles, for the origin 
of such a profile is still lacking, and there has been some consider- 
able debate over the exact functional form implied by the numerical 
studies. The fitting function most commonly applied has the gen- 
eral form 

pir^ccr-^^r. + ryo-^", (1) 
where 70 and 7„ are the power-law slopes of the central cusp 



and in the limit r — > 00, respecti vely, w hile is an appropriate 
scale radiu s. Dubinski & Carlbei^ jl99lh originally suggested the 
iHemq uis? ( 1990) profile, corresponding to yo = 1 and 7„ = 4, 
while[Navarro, Frenk & White (1996, 1997) argued for jq = 1 but 
7oc = 3 (the so-called 'NFW profile). Subsequently, s everal studies 



have argued for a somewhat steeper central cusp ("e.^. |Moor^^lJ 
1998lll999l;lGhigna et alll99i.l2000l;lFukushige & Makinoll99l 



20011) . with one recent suite of high-resolution simulations ap- 
pearing to imply 7o - 1.16 ± 0.14 (Diemand et al. 2004b). How- 
ever, even the highest-resolution simulation to date ( Diemand et all 
I2OO5I 70 = 1.2) can resolve the halo structure only to a fraction 
10"^ of its virial radius, and this for a single halo only. Much 
more common are numerical resolution limits several times larger 
than that, leaving room for the possibility that halo densities might 
become shallower than r~' at very small ('unobserved') radii 
<Tavlor & Navarrol200ll;|Power et all2003l ; lFukushige et all2004 
iHavashi et alj|2004) . and perhaps eveii tend to a finite density at 
r = Q with no cusp at all fStoehr et al.''2002"; "Navarro et al."2004|; 
iMerritt et al. 2005). Alternativ ely, there mig ht be no single, ' univer- 
sal' d ensity slope in this limit iNavarro et al . 2004; Fukushi ge et alJ 
|2004) . At the other extreme, there are very few hard constraints on 
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Figure 1. Plots of p|a■^^^ (left) and p/o"' (right) as functions of radius for ten simulated CDM haloes from lPiemand et all J2004jhl) . namely the six cluster- 
sized haloes called A9, B9, C9. D12, E9, and F9 and the four galaxy-sized haloes GO, Gl, G2. and G3. The radial coordinate in each halo has been scaled by 
the radius r2 at which dlnp/dlnr = -2 for that halo, and both ratios pI<tI„i and p/o"^ have been normalised by their respective values at r2. The lines in the 
upper panels are power-law fits to the scaled data from all ten haloes combined. The average negative slope and its uncertainty are indicated. In the bottom 
panels, we show the relative residuals [(p/o"^)-(p/cr^)fit]/(p/o"'). The rms residual is essentially the same (0.14-0.15) whether cTiot or o",. is involved, although 
in the latter case they are arguably distributed more uniformly about 0. 



any limiting value of the density slope as r ^ oo, as any halo is 
only well defined within a finite virial radius. 

In this paper, we examine the question of dark-matter halo 
structure from a slightly dilferent viewpoint, opting to derive p{r) 
from a simple dynamical Ansatz, rather than fitting a pre-set family 
of functions to simulated density profiles. Our starting point is still 
empirical, however, being based on another surprisingly uniform 
property of A^-body halo es. As was first noted and exploited by 
IXavlor & Navar ro (2001), the ratio of the density p(r) and the cube 
of the velocity dispersion a{r) is a single power-law in radius 

p(r)la-(rf oc r" (2) 

over the full numerically resolved radial range. IXavlor & Navarrol 
originally found a = 15/8 = 1.875, which interestingly is the value 
predic ted by the classic si milarity solution for spherical secondary 
infall iBertschineedl 1983) . They then used this constraint to solve 
the Jeans equation numerically (assuming an isotropic velocity dis- 
tribution) forp(r) and cr(r) separately. Although their result for p(r) 
is not of the form Q (it does tend to a shallow power law at the 
centre but steepens rapidl y outwards and falls to zero at a finite 
radius). Iravlor & Navarrol argued such a density distribution to be 
an adequate description of simulation data inside the virial radius. 
Unfortunately, there is no exact analytical expression for p{r) in the 
case a = 1.875, and this approach has not been used for fitting 
any data. Other studies have subsequently confirmed that p/cr^ is a 
power law in radius, but estimates of the exponent differ somewhat 
from the Taylor-Navarro value: a = 1.95 or 1.90 ± 0.05 according 
to|^sia, Tormen & Moscardini ( 2004) and Ascasibar et al. 1 2004), 
respectively. 

Our main aim is to investigate in more detail the structure 
and dynamics of spherical dark-matter haloes that follow a 'p- 
cr relation' of the basic type given in equation J2j, allowing at 



least initially for arbitrary values of a. Thus, in Figure \T\ we 
show this relation as defined by ten CDM haloes simulated by 
iDiemand. Moore & St adel ( 2004a b), the details of which were 
kindly provided to us by Jiirg Diemand. These include four galaxy- 
sized haloes and six cluster-sized haloes, with virial masses rang- 
ing from 10'- Mq to 10" Mq. All are dynamically relaxed. The left 
panels of Fig.Qexamine the ratio p/o"toi as a function of radius in 
these haloes, where the total one-dimensional velocity dispersion is 
crj^, = |(crj + o"a + o"^)- The right panels look at the quantity p/crj 
vs. r, where cr,. is the velocity dispersion in the radial direction only. 
Power-law fits to each of these profiles are drawn, and the residuals 
from the fits are shown in the bottom panels. These demonstrate 
that the ratio p/crj follows a power law in radius at least as closely 
as p/o-f^, does, although the fitted slopes differ slightly between the 
two cases. Fits for each of the ten haloes individually yield power- 
law slopes that can differ by ±4% from the average values in Fig.Q 
Whether this scatter is real or simply a reflection of numerical un- 
certainties is unclear, but it is certainly rather modest. 

To the extent that either of the p-cr relations illustrated in Fig.Q 
is 'universal', and insofar as dark-matter haloes are in equilibrium, 
imposing a dynamical constraint alon g the lines of equa tion to 
solve the spherical Jeans equation, as lTavlor & Navair ^ originally 
did, should lead directly to a 'universal' density profile. While we 
have no physical argument for the fundamental origin of the precise 
power-law behaviour in FigureQ(though clearly it must be related 
to the initial conditions and the formation via violent relaxation), 
it is much simpler to characterise than the density profile itself. 
Moreover, a Jeans-equation approach allows explicitly for a simul- 
taneous exploration of velocity anisotropy inside haloes — an issue 
which to date has been largely divorced from empirical descriptions 
of the halo density profiles. 

It is well known that the velocity distributions in dark-matter 
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haloes are not isotropic. We characterise velocity anisotropy using 
Binney's parameter 



Pir) = 1 - 



20-2 



(3) 



such that < yS < 1 corresponds to radial anisotropy and /3 < 
signifies a tangentially biased velocity distribution. It is typi- 
cally found that /3 ~ (isotropy) at the centres of haloes and 
gradually increases outwards (reflecting radial anisotr opy), reach- 
ing levels of j3 ~ 0.5 around the virial radius (e.g., IColm et all 
I2OOO'; 'Fukushise & Makino"200lV Indeed, it has been suggested 
(Cole & Lacev 1996; Carlbera et al. 1997) that there exists a 'uni- 
versal' anisotropy profile in dark-matter haloes . Very closely re- 
lated to such an idea is the recent claim by iHansen & Moord 
i2005h . that yS(r) depends roughly linearly on the logarithmic den- 
sity gradient 



7(f-) = -dlnp/dlnr 



with essentially the same, constant slope d/?/dy holding for a va- 
riety of end-products of violent relaxation processes (merger rem- 
nants, dark-matter haloes, collapse remnants). 

With these points particularly in mind, we base our analysis 
on the assumption that halo density and the radial component of 
velocity dispersion cr^ (rather than trtoi) are connected through a 
power-law relation of the general form 



Po 



ro 



(5) 



where ro is any convenient reference radius, po = p(ro), and CTro = 
cr,-(ro). Intuitively, as well as on the basis of FigureQ it seems most 
natural to expect the exponent 6 in this equation to be 6 = 3; but i t 
adds little complication to allow for the possibility iHansenl2004h 
that a slightly different value might provide a still more accurate 
description of simulated haloes. The choice of cr^ as the velocity 
dispersion to work with is arguably more natural than iXtot, since 
p(r) and 0", appear separately in the spherical Jeans equation but 
cr^jjj = 0-2(1 - 2/J/3). Technically, this choice ultimately allows for 
a more tractable inclusion of velocity anisotropy; empirically, it is 
obviously well justified by FigureQ 

We begin in Section |2| with an investigation into which 
spherical density profiles satisfy the Jeans equation and obey the 
Ansatz under the restrictions of velocity isotropy (J3 = and 
o-,ot = o",-) and a fixed 6 = 3. The problem is t hen identical to the 
one first considered bv lTavIor & Navarrol i200 1 ). and our approach 
is rooted in theirs, but we a l so dra w on some aspects of the consid- 
erations by I Williams et alj i2004) . However, unlike these authors, 
we explore the full solution space of the problem. We find that only 
very few of the many possible solutions correspond to realistic den- 
sity models for simulated dark-matter haloes, and in fact only one 
solution, which occurs for a 'critical' value of a, is of practical im- 
portance. 

We then proceed in Section to consider the more realis- 
tic case of anisotropic velocity distributions and allow for gen- 
eral values of 6 in equation j5}- We show that in the case of an 
anisotropy parameter p that depends linearly on 7 (including con- 
stant anisotropy as a special case), and for any e, the solutions of 
the Jeans equation under our adopted constraint are exact analogues 
of those in the /? = 0, e = 3 case. In particular, for each pair 
(e, j8(r = 0)) only one physical solution of practical relevance ex- 
ists, which again occurs at a 'critical' a value. These solutions have 
fully analytical density, mass, and velocity-dispersion profiles with 
power-law asymptotes at small and large radii. 



In Section |4| we compare our analytical profiles to the 'ob- 
served' density, velo city dispersion, and an isotropy profiles in 
haloes simulated by Diemand et al. ( 2004a b, 2005), finding good 
agreement in general. Finally, Section|5|discusses our findings and 
summarises the paper. 



2 THE ISOTROPIC CASE 

Our underlying assumption is that some version of the general 
p-cr relation in equation Js} holds for dark-matter haloes. Before 
allowing for this level of generality, however, there is much in- 
sight to be gained from beginning with a more specialised case, 
in w hich the velocity d istribut ion is isotropic and e = 3. Then, 
as in Tavlor & Navarrol f200lh . lwilliams et alj <2004 . and lHansenI 
(i2004) . we have 



(4) A = ilL I J_ 



,3 



(6) 



The Jeans equation for a spherical, self-gravitating collisionless 
system with isotropic velocity distribution is 

d(pcrj) ^ ^ GM{r) 



dr 



-P(r) ■ 



(7) 



with M{r) = An J^' ip-p{u) du. Following iTavlor & Navarrol we 
solve equation ^6} for a,-, insert it in equation 0, and differenti- 
ate again to obtain 

[x^ d 



1 d 

x^y dx 



y dx 



Here, x = r/ro and y = p/po are dimensionless variables, and 
K = AnGporl/a-lg 



(8) 



(9) 



is a dimensionle ss measure of the velocity dispersion scale. 
iTavlor & Navarrol studied the solutions of equation Jsj for the par- 
ticular value tt = 15/8 by numerical integration. It proves useful, 
however, to first re-write the problem in terms of the (negative) 
logarithmic density slope yir), as defined in ^4). Equation Jsj then 
reads 

/ - f (7 - ja) (r - n) = Ikx--^"'V\ (10) 

where a prime denotes differentiation with respect to In x and 

ya = \a, Jt = a+l- (11) 



The objective of this Section is to investigate the solution 
space of equation <10> . First, note that its r.h.s. becomes constant 
for y = X"''' , with 

71=6-20'. (12) 

In fact, as was already noted bv ITavlor & Navarrcl this singular 
density profile is a solution to equation <I0> and, for a = 2 cor- 
responds to the well-known singular isothermal sphere. For 3/2 < 
a < 5/2 (which covers our regime of interest; see FigureQabove), 

Ja <7i < Jb (13) 

and thus, following ITavlor & Navarrol we choose to identify the 
reference radius ro in equation jSJ as that at which the (negative) 
density slope equals 71, i.e., 71 = y{x = 1). In this way, ro is well- 
defined for all realistic solutions' . Moreover, it is then obvious from 

' As we shall see, for some solutions y{r) = 71 at more than just one radius. 
In this case, we pick that radius which in addition maximises y' . 
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a = 15/8 




a = 35/ 10 
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Figure 2. Flow diagrams of equation 1141 . describing isotropic density profiles with p/cr^ oc r"", in the (y, y') phase space (the length of the arrows has been 
normalised for ease of reading) for three dilferent values of a. The flow vanishes at the three fixed points (dots), where y' = and y is a root of the light-hand 
side of equation 1141 . Some solutions to equation 1141 are given as curves. Physical solutions to the original equation 1101 must lie above the (dashed) parabola, 
which is defined in equation 1151 . Solutions with y' < or y < at any point (red and magenta curves) are unrealistic. There is exactly one realistic and 
physical solution for every a ^ 35/18 (blue curve) and none for larger a. At small radii this solution's density profile approaches the power-law p oc r"^"'^, 
while at large radii the density is truncated at a finite radius for a < 35/18, and approaches p oc r^"^^!^ = ^^31/9 tor a = 35/18. 



this equation that the constant k is effectively a measure of y' at 
X = 1. 

Let us now consider equation <10> in some detail, in particular 
the possi ble behaviour of any solutions at very small and large radii 
(see also lwilliams et all) . It is straightforward to verify that in the 
limit r — > only three asymptotes are possible: either y — » 
(both sides of the equation vanish, so long as ff < 5/2); or y — > 
7i (both sides approach a constant); or the density has a central 
hole with y = inside a non-zero radius (note that y ^ y;, is not 
possible, since the r.h.s. would then diverge for a > 3/2)} Of these 
three, the latter is unrealistic, whereas y — » yi is inconsistent with 
simulated haloes, which have shallower inner density slopes. This 
leaves y — > y^ as the only interesting option. We thus expect that, 
for any given value of a, there exists one solution with y — > y„ as 
r — > 0. The solution with this limiting behaviour is associated with 
a unique value of k, which depends on a through equation <10> . In 
fact, this very case corresponds to the 'critical k' solution developed 
bv lTavlor & Navarrd for their specified a = 15/8. 

In the limit r — » oo, also only three possible asymptotes exist: 
either y — > y;, (both sides of the equation vanish, if a > 3/2), or 
y — > yi (both sides approach a constant again), or the density has 
an outer truncation with v = beyond a finite radius (note that 
y — » y„ is not possible, as the r.h.s. then diverges if o- < 5/2). All 
of these three are physically meaningful (apart from the fact that 
for y — > yi the mass forma lly diverges, similarly to the isothermal 
sphere). fTavlor & Navarrol ' s 'critical k' solution for a = 15/8 is 
one with an outer truncation. 

In order to gain more insight into the solution space of equa- 
tion <10t . we now follow IWilliams et alj and differentiate it yet 
again (with respect to In x) to obtain 

r" - r' (r - |[2y<, + 2yi, - y,]) = f (y - yj (y - y*) (y - yi), (14) 



2 iHansenl 12001 . in his analysis of the same problem, mentions the asymp- 
tote y — > yi in the Hmit r — ► 0, but misses the important possibility y — ► yo- 



equivalent to their equation (2.2). In Figure |2| we plot flow dia- 
grams of the (y, y') phase space for three representative values of 
a. The dashed parabola in each panel corresponds to 

r' = Tmin =j(y- 7a) (j - 7b)- (15) 

While this actually is a solution to jl4t . as already discovered by 
Williams et al., it is not a valid solution to the original, physical 
problem. This is because the constant /c > by definition, so equa- 
tion not requires that y' > y'^^^^^ in order for its r.h.s. to be non- 
negative. In fact, the equality y' = y'^-^^ is only allowed in either 
of the limits y — > y„ as r — > or y — > yj as r — > cx) (see the 
discussion above). Similarly, any solutions below the parabola <15t 
require k < and are unphysical. By contrast, all solutions above 
the parabola <15> correspond to potentially viable solutions of the 
original equation llUt . 

In all panels of Figure |2| the flow vanishes at the three fixed 
points (y, y') = (ya,0), (yi,0), and (y*, 0), which are plotted as 
blue dots. The central fixed point is stable and corresponds to the 
singular solution y = x"''' discussed above. Solutions with y" > 
(arrows pointing upwards) at any y > yj approach, in the limit 
y — > CM, the upper right branch of the parabola y' = y[^^^- This 
corresponds to the case, discussed above, of an outer truncation to 
the density. More precisely, integrating equation <15> twice yields 
the asymptotic density profile 

p _ ,-r.(^2(r*-r»)/3 _ ^2(y,-,„)/3)3/2 ^jg^ 

as r approaches some finite radius r,. For r > r,, p = 0. Con- 
versely, solutions with y" < (arrows pointing downwards) at any 
y < yi approach, in the limit y — > -oo, the upper left branch of 
the parabola. This corresponds to the case, mentioned above, of an 
inner truncation, with the asymptotic density profile 

p ^ r-^^irfy-y"*'^ - ^2(r»-)'*)/3)3/2 

as r approaches some radius r;, > 0. For r < f;,, p = 0. Thus 
solutions for which y — » cx) at large radii have a finite outer trunca- 
tion radius, and solutions for which y — > -00 at small radii have a 
finite-sized inner hole. These latter solutions are clearly unrealistic 
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and also unphysical (the isotropic distribution function must be- 
come negative to account for the hole). Furthermore, solutions that 
ever visit 7' < are hardly realistic (their density profiles become 
shallower towards larger radii, at least over some radial range). 

Apart from this — as Williams et al. (2004) have also noted — 
the generic behaviour of the solutions to equation <14> depends on 
whether a is greater than, less than, or equal to a critical value o-crii 
for which the inner and outer fixed points 7,, and yt are equidis- 
tant from the central 71. Referring to equations II H and 1121 . 
ri = \{ya + Jb) requires a^n = 35/18 = 1.94. 

The behaviour of the solutions for a > cfcrit is exemplified 
in the right panel of Fig. |2| The solutions are split into two fam- 
ilies by the one (upper magenta) which ends at the fixed point at 
7 = 7ft, so that p oc r"^* as r — » cx). Solutions lying above this curve 
correspond to density profiles with a central hole and an outer trun- 
cation; those below usually also have a central hole but at larger 
radii perform damped 'oscillations' about 7 = 71, eventually ap- 
proaching p oc r"''! as r — » 00. All of these solutions are unrealistic, 
since they possess inner density holes. The only exception is a lim- 
iting solution (lower magenta) which starts from the fixed point at 
y = fa (i.e., p oc r"''° as r — > 0) and slowly approaches 7 = 71 as 
r ^ 00. However, this solution is still not a viable description of 
dark-matter haloes, since its density profile is too shallow at large 
radii. 

The behaviour of the solutions when a < a^ni is illustrated in 
the left panel of Fig.|2|for a = 15/8, the case previously studied by 
[javlor & Navarro. The situation is in some sense a mirrored version 
of that for a > ttcrit- The solution space is again divided into two 
families, now by the solution (blue) starting from the left fixed point 
at 7 = 7a at r = 0. Solutions above this one again correspond to 
density profiles with a central hole and an outer truncation; those 
below it start from an 'oscillation' about the power-law p oc r"''' 
in the limit r 0, eventually steepening outwards and generally 
being truncated at a finite large radius. The limit of this family is 
the solution (magenta) which instead of an outer truncation has a 
power-law fall-off with p oc r""' as r — » 00. 

For any a < ci-crit, the (blue) solution separating the two fam- 
ilies just described is a potentially realistic model for dark-matter 
haloes, starting as it does from a shallow power-law cusp p oc 
at r — > (with y^ = 2a/5 < 7/9 for a < 35/18) and steepening 
outwa rds to reach p = a t a finite radius. This is the model identi- 
fied bv lTavlor & Navarrol as their 'critical k' solution in the specific 
case a - 15/8. Given this a, numerical integration of equation <14t 
outwards from 7 = 7^ = 3/4 at x = yields k = 2.61 A, to be com- 
pared with the value k = 2.619, found by Tavlor & Navarro through 
trial-and-error integration of equation jsj starting from x = I. 

The middle panel of Figure|5|depicts the special situation a = 
tj-crit, for which the flow in (7,7') phase-space is symmetric with 
respect to the transformation r — > r"' and 7 — > 271 - 7 (because 

7 a + Jb ■■ 

integral. 



271 by definition of ffcrit) In this case, there exists a first 



K 



(/ + |(r - r«)(r - r*)) (r' - f (r - roXr - i'b)f 



(18) 



which is conserved by any solution of equation I14> with a = a^, 



In particular, K 
also for 



for the solution 7' = 7j„i„(7) (see eq. 1151 ) and 



/ = -|(7-r<.)(r-7*) (19) 

which is plotted as the blue curve in the middle panel of Figure|2| 
It again divides the solution space into two, although now solutions 
above the curve )19t always have both an inner hole and an outer 
truncation to the density profile, while solutions below the curve 



undergo undamped 'oscillations' about p oc r"''', never settling to 
an asymptote in either limit r — > or r — > 00. 

The solution <19> for a = a^^n = 35/ 1 8 is particularly appeal- 
ing since it starts from a shallow power law 7 — > y„ = 1/9 in the 
limit r — > and tends to the steeper 7— »7(, = 31/9asr— >cx). The 
combination of these features is reminiscent of dark-matter haloes. 
This single solution with a = ci-crit is the only one, of all the solu- 
tions for any a, whose density has both inner and outer power-law 
asymptotes and monotonically increasing 7. It is also the only one 
we have found which is simple enough that almost all its physical 
properties can be developed analytically: equation <19> is easily in- 
tegrated to give y(x) and subsequently y(x) oc p(r) as simple func- 
tions, which then allows us to evaluate o-,.(r), the enclosed mass 
profile M(r), and the circular velocity Vc(r). 

To aid in obtaining these basic results, we first substitute equa- 
tion 1191 into equation jlOt and evaluate the result at x = v = 1 
(where 7 = 71 by definition) to find 



^= fi(ri-r«)(r/.-7i) = f 



(for a = ffcrii)- With this in hand, we obtain 
5M„ 



P(r) 
y(r) 



o-j(r) 



-1/9 



9^rl 



(1 + x^'^y 



l + x^/s ' 
9 GM,„, 



40 ro 



l+X^f 



M(r) = M,„ 



1 + x^/' 



V^(r) 



ro 



^4/9 



l+x^/' ' 



(20a) 

(20b) 

(20c) 
(20d) 
(20e) 
(20f) 



where x = r/r^ as usual. We have replaced po with the total mass 
Mot = (576;r/5) po'"o; and the value of k in equation J20at has been 
used with the basic definition J5) to eliminate the normalisation (T,-o 
from equation <20dt . Note that, because this solution has a finite 
total mass, both tr^(r) and V^(r) fall off as r"' in the limit r — > 00. At 
the same time, both quantities vanish at r = 0, and thus each profile 
peaks at a finite radius. This happens at x = (7/9)''"* for cr^(r), 
and at X = (11/9)"'* for V^(r). Finally, the gravitational potential 
follows from integrating dO/dr = GM(r)/r^: 



<^ir)-l^ 
4 ro 



B 



(20g) 



whe re B„(p, q) = C I '' '(1-0* 'df is the incomplete beta function 
(e.g. lPress et alJl992l §6 .4). For an alternative form and asymptotic 
limits of C), see equation j40i> and the following text. 

Apart from its convenient and unique analytical properties, 
this solution of the Jeans equation is additionally of interest because 
it corresponds to a p-cr relation of the type in equation ^6} with 
o: = ffcrit = 1.94, remarkably close to the exponent actually found 
for the simulated dark-matter haloes shown in Fig.Qabove. An ob- 
vious caveat is that our development in this Section has assumed 
an isotropic velocity distribution, which is known to be incorrect. 
Nevertheless, having characterised the solution space of equation 
^H) or <10> for this restricted case, it turns out to be straightforward 
to allow for realistic velocity anisotropies and at the same time in- 
vestigate other values of e in the full p-cr relation of equation js}. 
As we will now show, the specialised isotropic solution of equa- 
tions <20> is in fact just one of a larger, more general family of 
analytical solutions to the Jeans equation. 
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3 THE ANISOTROPIC CASE 

We now return to the more general form of our basic assump- 
tion namely 



and seek to solve the spherical Jeans equation in the general form 
d .mr) 2 , , GM{r) 

- pa,. + —pa^= -p{r) (21) 

with anisotropy parameter /? as defined in equation J3j. By the same 
method as in the last Section, we then find the generalisation of 
equation <10t to be 



r - 



26 

2 + e 



3')-^(r--|)(r-5|^/^- 



20- 
2 + e 



2 + e 



(22) 



Differentiating again to obtain the equivalent of equation <14> leads 
to a number of nonlinear terms involving up to the second deriva- 
tive j8", including cross-terms of the type yP' and ySy'. These terms 
cancel exactly, however, if and only if p depends linearly on y. That 
is, the structure of the Jeans equation itself naturally suggests that 
we stipulate the relationship 



= + b(y - y„) , 



(23) 



for b and constants. The definition of y^ follows from substitut- 
ing the expression into equation <22> . yielding 



2 6 

r' - - (r - 7a) (r - 7b) = ^r- t-j ' 

6 2 + e - 2eb 

if 



7o 



2a 2e 

+ fia- 



(24) 



(25) 



2 + £ 2 + £ 

Here x and y are defined as in Section|2| while k is given in equa- 
tion J3 and 

6 



7h 



(26) 



Equation <24> is thus the generalisation of equation <10> for arbi- 
trary e and p linearly dependent on y (for 6 = 3 and Po = b = 0, the 
two equations are identical). Note that b < (2 + e)l2e is required to 
keep the r.h.s of equation <24t non-negative. 

Another differentiation of equation 1241 with respect to Inx 
yields the generalisation of equation <14> 



6 - e , / 27„ -t- 274 - (6 - 2)7i 

r r 7 

6 \ 6-6 

2(^-2) , 

— 5 — (7 - 7a)iy - rfc)(r - n) 



with 



71 



2(6 -g) 

6-2 



(27) 



(28) 



Again, the singular function p oc r '''' is a solution, also identified 
bv iHansea C2004) . As long as 



e> 2 



(29) 



to avoid fundamental changes in sign, the only differences between 
equations <27> and <14> are in numerical values of constants in- 
cluding 7„, yt, and y\. The topology of the (7,7') phase-space for 
general 6 > 2 with anisotropy parametrised as in <23> is thus the 
same as for the isotropic case with 6 = 3, and mathematically the 



solution spaces of the two problems are isomorphic. Hence, most 
of the discussion around Figure|2|carries over here. In particular, 
for any given 6 > 2 and /Jq < 1 , the fixed points of equation I24> at 
7fl and 74 satisfy 7^ < yt, and they bracket the third fixed point at 
71 as long as 



6 2-1-6 

- <a < 

2 2 



(30) 



As before, there then exists a critical value of a for which the three 
fixed points are equally spaced: with the requirement 71 = \{ya + 
yi,), equations <25t . <26> . and <28> give 



(10-6)(2 + 6) 2(6-2) 



2(6 + 6) 



6 + 6 



■Po 



(31) 



which notably does not depend on the slope b in the linear p-y 
relation of equation J23> . With 6 = 3 and po = 0, we have = 
35/18 as in the previous section. 

The generic characteristics of the various of solutions to equa- 
tion J27t are again determined by whether a < a^-^n, a = a^^-n, or 
a > Q-crii (see Figure|2j. For any constant anisotropy {fio <, i and 
& = in eq. 1231 ). the division between physical and unphysical 
solutions in all three cases is exactly analogous to the isotropic 
specialisation of Section |2| When a gradient in p(r) is allowed, 
however, the number of physical solutions to the problem becomes 
smaller, because p is bounded above by 1 at all radii. In particular, 
for the most relevant case ofb > Q (corresponding, for realistic halo 
models, to an outwards increasing radial anisotropy), 7 must also be 
bounded above, and solutions with an outer truncation to the den- 
sity profile (7 — > 00) are no longer viable. As a consequence, the 
only physically possible solutions with a shallow density cusp in 
the centre and monotonically increasing slope 7 are the analytical 
solutions that occur only for a = o-crii. 

Setting a = a^f^ in equation )27t reduces it to 



„ 6-6 , 2(6-2) 
7 7(7-71)= 5 — (7-7a)(7-7f.)(7-7i) (32) 



with 



10-6 2(6 + 2) 

7a(Q'cril) = + Pa, 

6-1-6 6 + 6 

10 + 36 4 

7l(Q'crit) = — + Po, 

6+6 6+6 

10 + 76 2(6-2) 
7i(a'ciit) = 2 P^- 

6 + 6 6 + 6 



(33a) 
(33b) 
(33c) 



Generalising equation <18> . any solution to equation <32t conserves 
a first integral 



^ = (r' + #(7 - 7«)(7 - 7*)) (7' - f (7 - 7«)(7 - 7*))"*""- (34) 



Thus, for A" = in particular two simple solutions exist. One is 
y' ~ y'mm = f (7 - 7o)(7 - 7b), which defines the parabola in (7,7') 
phase-space within which all physical solutions must lie (cf. Figure 
|2j. The other is 



f- 2 

7 = - -7. — (7 - 7a)(7 - 7b), 
2e 



(35) 
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which is easily integrated to give an analytic expression for j(x) oc 
p{r)? Before writing down this and other quantities, it proves useful 
to introduce the auxiliary parameter 



'-{jb-Ja) = 2 



(e-2)(2-/?o) 



2e 6 + e 

the meaning of which becomes clear below. Then 

^4 + 3;7-2j8o 



4 - 7/ - 2^0 
so we have from equation j31> that 
477 

ttcril =T} + 2- — , 

4 - ?7 - 2/3o 

and from equations <33> . 

To = ra(acrit) = 1 - '7/2 +/3o, 

71 (a,,,,) = 2 + 77/4+/3o/2, 
Too = r6(»ciit) = 3 + 

Integrating equation <35> . we then find for the density 



(36a) 
(36b) 

(36c) 

(36d) 
(36e) 
(36f) 

(37) 



Thus, the parameter 77 governs the speed of the transition between 
the power-law asymptotes p oc r^f^ at small radii and p oc r"''~ 
at large radii.** In addition, the constant /3o in equation <23> takes 
on physical meaning as the velocity anisotropy at the centre of the 
density distribution. 

Astonishingly, the gradient b = AjijAy in equation <23> does 
not appear in equation <27> or any subsequent relations, implying 
that the influence of anisotropy on the density profile is entirely 
determined by the situation in the centre, if [} depends linearly on 
y as we have assumed. For any fixed e > 2, the main effect of 
a radially biased velocity ellipsoid at the centre (0 < ySo ^ 1) is 
to steepen the inner power law 70 relative to its isotropic value, 
and make 700 smaller, i.e., the outer density profile shallower. The 
reverse holds for a tangential anisotropy, pQ < 0. In order to keep 
7o > then (so the density does not decrease towards r — > 0), we 
require 



poXe- 10)/(4 + 2e), 



(38) 



which excludes very strong tangential biases. Isotropic models are 
allowed only for e < 10, a limit which is not relevant to dark-matter 
haloes, but corresponds to the classic Ipiummed il91lh sphere. On 
the other hand, the physical requirement pQ < I implies 70 < 2 
and 7c<, > 3 for any e > 2. As a result, all the critical-o- density 
profiles in equation <37t have a finite total mass, and fully analytical 
cr^(r) and M(r) profiles which are of the same basic form as in 
equations <20> . 



^ Note that in the Hmit e = 2, equation 1351 is simply y' = 0, while equation 
1311 gives ffciit = 2 for any /?o so the the right-hand side of equation 1241 
becomes a non-negative constant. Therefore any of a continuum of singular 
solutions, y = x^''' for y, a constant in the interval ya < < ji,, is a va lid 
critical-Q- solution for e = 2. This is the situation that iHansenl |2004) is 
strictly relevant to, although pure power-law density profiles such as these 
are not applicable to simulations of dark-matter haloes. 
"* These sol utions are me mb ers of the much broader class of 'ojSy' models 
discussed bv lZhacllT993) . In lZhaoT s notation (which is completely dift'erent 
from ours), the profiles of equation 1371 have (a.yS, y) = ('/"', yoo, 70) = 
(;;-',3 + 77,1- 77/2 +/?o). 



To obtain these profiles in detail, we first express the linear 
relation between p and 7 in terms of p^ = Vimr^„ p(r) instead of b, 
namely, 

Poo-Po 



P(r)=Po + ' 



■ (y(r) - Jo)- 



Too -ro 

Then it is straightforward to show that for 

1 



P(r) 
y{r) 
P(r) 



-(4-h77-2/3o)(4-h77-2/3„), 
o 

4-H77-2/3oM,ot ^^__^„^^ + ^y(y^-y«)h,_ 



7o + 7°°-^'' 
\+ x'^ ' 

1+ x'l ' 

1 GM,, 



y \ 1 



4 + 1]- 2y8„ ro \ I -I- x'l 
M(r) = M, 



1 -I- x'l 



x^ 



I -I- xi 



(39) 

(40a) 
(40b) 

(40c) 
(40d) 

(40e) 
(40f) 
(40g) 



where M,o, is related to po by the condition p(7- = ro) = po and 77, 70, 
and 7„ are, of course, given in terms of e and Po by equations j36> 
and satisfy (7„ - yo)/77 = 26/(e - 2). The total one-dimensional 
velocity-dispersion profile in these models is also analytical, being 
given simply by o"^o((7') = (T^(r){l - 2p(r)/3). Note that for e = 3 
and y8o = ySoo = (the specialised case considered in Section 
77 = 4/9 and we recover all of equations )20t . As in that case, the 
velocity-dispersion and circular-velocity profiles show peaks in this 
more general situation: crj has its maximum at x = (1 - 77/2 -/Jq)'''' 
at X = (1 + 77/2 -j8o)''''. Finally, the gravitational potential 



and V^, 
is 



0(7-) = - 



r]ro 
rjro 



TTW \1 'I 27 

V';' 1 27 \ ri 2' r;7J 



(40h) 



(40i) 



where again Bu{p,q) = ^ ?''"'( 1 - f)' 'df is the incomplete beta 
function and B(p,q) = Bi(p,q) = r(p)r(^)/r(p + i5') the (complete) 
beta function. In the limit of large radii O — > -GMiotT""' , while for 
small radii 

0(.) . [iB{i, . i) - . (40j) 



'■0 



It is worth noting that at small radii p oc r '"^ and crj oc . Thus 



the pressure pcr^ 



, which diverges for radial anisotropies at 



the centre (Po > 0) but approaches a constant for central isotropy 
and vanishes for tangentially biased central velocity distributions. 

Aside from its pleasing — and somewhat surprising — 
simplicity in the face of a nontrivial radial variation of P(r), the 
family of models defined by equations )40t is further interesting 
because the linear 'P-y relation' in equation j23> or <39> is pre- 
cisely what iHansen & Moord i2005h have suggested is a generic 
result of collisionless collapses, mergers, and relaxation processes 
(see their Figure 2). 

Moreover, it is beneficial that for any e these anisotropic mod- 
els have the same critical value of the exponent a in our p-cr relation 
jsj as do the fully isotropic models, just so long as isotropy holds 
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at the centre alone (/3o = 0). Simulated dark-matter haloes indeed 
tend to be roughly isotropic at their centres and radially anisotropic 
in their outer parts. Thus, as was discussed at the end of Section|2| 
it is again striking that the data shown in the right-hand panels of 
Figure Q exhibit a scaling p/cr^ ~ r^" with a ^ 1.92 very close to 
the expected critical value (35/18 = 1.94) for 6 = 3 and j8o = 0. 

With these points in mind, we now go on to a detailed 
fitting of our models in equations <40> to the density and 
velocity-dispersion profi les of dark-matter haloes simulated by 
iDiemand et^j2004Jhft . 



4 COMPARISON WITH SIMULATED HALOES 

In order to compare the anisotropic models in equation <40t against 
numerical dark-m atter haloes, we again make use of the simula- 
tions published bv lPiemand et alJ <2004a b), which we referred to 
in Section Q (Fig. 0. To repeat, these include four galaxy-sized 
haloes (virial masses 1 - 2x 10'^ Mq) and six clusters (virial masses 
2.4 X lO"* Mo - 1.3 X 10'' Mo). All were evolved to redshift z = 
except for two cluste rs which w ere run further ahead in time to 
complete mergers; see lPiemand et aO for full details. 

In addition to density profiles, the data on these haloes include 
the enclosed mass profiles M(r) and the radial and tangential com- 
ponents of velocity dispersion, a^^(r) and cr^(r) = o"g(r) -I- cr^^{r). 
From the run of p(r) in a series of spherical shells we have 
estimated the local density gradient as 7(r,) = -[log p(r,+i) - 
logp(r,_i)]/[log - log From the velocity dispersions, we 
have calculated ;S(r,) = 1 - cr^(r,)/2cr^(r,) and the total one- 
dimensional dispersion, arlaSn) = |[o"r(''i) + o"?('"i)]- 

Five parameters define the models in equations <40> : the nor- 
malisation constant Mtot and the scale radius ro; the parameters e 
and Po, which fix rj, yo, and y„ (and thus the shapes of all profiles); 
and finally the velocity anisotropy at r — > oo. Together these 
must suffice to describe the separate p(r), a--^(r), and p(r) profiles 
for a dark-matter halo. While it is of course possible to fit each of 
the ten simulated haloes individually, we are more interested here 
in the question of whether some 'universal' parameter values might 

apply- 
Thus, we start by assuming that all haloes follow the same p-cr 
relation, p/cr,': oc r"""'', with a single value of e. That this is likely so 
is already suggested by Figure □ We further assume that the central 
velocity anisotropy po is the same for all haloes. This is more of a 
debatable contention — even though most simulations show similar 
(low) levels of anisotropy in their innermost resolved regions, there 
is no clear evidence that po always tends to a single value. Never- 
theless, if e and Po both are 'universal' then the shape of the density 
profile must be too in these models, and so this possibility is worth 
examining. 

We proceed by defining a grid of (e,jSo) values. For each pair 
in turn, we compute t/, yo, and yon from equations <36> . which al- 
lows for the calculation of dimensionless model profiles Pmaiir) and 
<Tjj^^j(r). For each of the ten haloes, we then find values for Mjj,, 
r'■^^\ andpf£ to minimise the sum of absolute deviations, 

A;(e,/3o) = }A log — r-- + log r— \ . (41) 

jr^ { pin) o-r(n) j 

which is more robust against outliers than the standard statistic. 



Here Nj is the number of data points in the jlh halo; typically, Nj = 
21. The total deviation that is minimised for each (e,Po) is thus 

A,o,(e,y8o) = Yj ^(f'/^o) • (42) 

Ultimately, we find the set {e,^o; M,'^*, r'g^\p'£\ 7 = 1, . . . , lOj for 
which Atot is the minimum over our original grid. Strictly speaking, 
this occurs at e = 3.2 and Po = -0.05, but the minimum is rather 
shallow and the best fit with e = 3 exactly (which requires Po = 
-0.10) is not significantly worse. The details of this latter fit are 
shown in Figure^ 

The upper panels of Figure show the density and (radial) 
velocity-dispersion data that were used to constrain the model pa- 
rameters. To emphasise the shapes of the distributions, the radial 
coordinate in each halo j has been normalised by the fitted value 
of rg \ and the densities and velocity dispersions have been nor- 
malised by their fitted values at r'^\ The best-fitting dimension- 
less profiles, from equations <40b> and <40et . are drawn as the bold 
curves. The rms relative deviations from these curves, (p -pmod)/p 
and (cr; - crj^^,^)/o'l, are both of order 1 1%-12% for all ten haloes 
combined. 

The left-middle panel of Figure|^shows the negative logarith- 
mic density gradient y(r) as a function of r/r^^ in the simulated 
haloes against the model curve given by equation )40c> for 6 = 3 
and Po = -0.10. Given these parameters, equations <36t imply that 
the power-law slopes at r — > and r — > 00 are yo = 2/3 and 
y„ = 52/15 = 3.46, respectively. The transition from the inner 
to the outer power law is rather gradual, with = 7/15 = 0.66. 
Recall that ro is defined as the radius at which 7(^0) = 71, with 
7i = 5(70 + 7oo) = 31/15 in this case. 

The right-middle panel then shows the 'observed' anisotropy 
parameter p(r) vs. the scaled radius x = r/r^K The bold curve 
traces the model relation j40d> for po = -0.10 and jS™ = 0.67, 
which is the average of the ten different p£ values obtained by fit- 
ting the (tJ profile of each halo separately. The dashed curves have 
the same Po = -0.10 but p^c = 0.53 and p^c = 0.84, correspond- 
ing to the minimum and maximum of the fits to the ten haloes. 
Evidently, these curves together account for much of the observed 
scatter in p{r). Note that the average /?„ value implies a slope for 
the p-y relation 03 of (j8„ -po)/iy„ -70) ^ 0-28. This IS compara- 
ble to the slope 4e/d7 » 0. 19 inferred bv lHansen & Moore! 1200^ 
from their investigations of a completely different set of simul ated 
haloes. However, it should also be noted that iHansen & Moord 
found a relatively tight correlation between p and y (see their Fig- 
ure 2). Comparison of the data points in the two middle panels of 
our Figure shows somewhat more scatter in any empirical p-y 
relation for the haloes we are working with. This is particularly ev- 
ident at relatively large radii, r/ro > 2. In fact, it is not clear that 
a 'universal' slope in a linear p-y relationship (i.e., a unique value 
of Poo ) can describe all of these data if the density profile is strictly 
'universal' (i.e. if there truly is a single value for Po as well as e). 

The bottom two panels of Figure |3| complete the comparison 
with our models. The left panel shows the re-normalised circular- 
velocity profile, ^JM{r)rolM{ro)r vs. r/ro, which is essentially 
equivalent to the top-left panel showing the model and observed 
density profiles. The right-hand panel shows the normalised to- 
tal one-dimensional velocity dispersion profile, which is obtained 
from the two panels above it: crl„f{r) = crj(r)(l -2/?(r)/3). Note that 
the scatter of cr^^j, about our model is some 25% smaller than that of 
crj: the scatter away from the model curve in the upper-right panel 
is compensated by the scatter of /3(r) in the middle-right panel. 
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Figure 3. Comparison of our analytical, ' critical-cc' models in equation 1401 with the same ten simulated CDM haloes from lPiemand et all l2004all]|) that have 
already been employed in Fig.Q The exponent e in the p-o" relation js) was kept fixed at e = 3, while a global fSg was fit to the density and velocity dispersion 
profiles of all haloes simultaneously. For each halo, individual values for the scale radius ro, normalisation constant M^at, and asymptotic anisotropy f!oo at 
r —> CO have been fitted (see text for more details). The solid curves show the dimensionless model profiles given the global e and /Jq. Squares represent the 
data from the simulated haloes, with the various profiles for each halo scaled by their best-fit values at the fitted ro for that particular halo. For the anisotropy 
profile, the solid curve represents a mean model, in the sense of a mean value for /?„, while the dashed curves indicate the range of best-fit values for/3oo. 



Overall, it is surprising how well our simple model is able to 
reproduce the main features of the spatial structure and the kine- 
matics of these simulated dark-matter haloes. However, it remains 
to be checked that the p-a relation in the 'observed' haloes is con- 
sistent with that required by the analytical models we have fit. 
Specifically, by using equations )40t we have assumed that the 
value of a is the critical one given by equation <31> . For 6 = 3 
and j8o = -0.10 as in Figure|3l this is o-cdt = 59/30 = 1 .96. Fig- 
ure|^plots the p/crj data points from the ten lDiemand et a l. haloes 
(normalised by our fitted values of p^ ' and cr* ^ for each halo indi- 
vidually) against the scaled radius r/rj,^'. The dashed line has slope 
-ttcrii and gives a reasonable description of the data. In fact, direct 
linear regression (with a 3-it clipping applied) yields a best-fit slope 
of -1.933, only^l.5% different from the expected value.^ This is 



^ The fitted slope here differs slightly from that in the right-hand panel of 
Fig. n because now we have scaled to the fitted radius r^^ in each halo, 
rather than to model-independent, but cruder, estimates of ri . 



drawn as the solid line in Fig. |4] The relative residuals from the 
critical-cf model line, [(pjcr].) - (p/o-l)„^oi]/(p/crl), are shown in 
the lower panel. Evidently, the largest deviations from this power 
law occur at r/rg > 2, which is where the crj(r) and p(r) profiles 
scatter most in Figure|3| 

As a last test, we wish to compare our model against the es- 
timated densi ty profile only of an e xtremely high-resolution halo 
simulated_bv jDi£mand et aP i2005h . This is again a cluster-sized 
halo. As lDieinan^^n describe in detail, this particular halo was 
defined by first performing a simulation with spatial resolution of 
10"^ in units of the virial radius, evolving the run to a high redshift. 
The most central part of the density profile was then scaled to match 
onto the outer parts of a lower-resolution halo previously evolved 
to z = 0. Because of this estimation procedure, we do not have the 
full velocity-dispersion and anisotropy profiles for this cluster. 

Even in the absence of kinematical data, we can fit for all 
of e, fio, M,o,, and using the density profile alone. In prac- 
tice, we set e = 3 and only fit for j6o and the two normalisa- 
tion factors by minimising the sum of absolute deviations A = 
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Figure 4. Top panel: plot of p/crj against radius for the ten simulated CDM 
haloes already employed in Figs.Q]and|3] In fact, this plot is very similar to 
the right panel of Fig.0 the only distinction being that we use our best-fit 
values for individual scale radii rjj'' to normalise the data in each halo. The 
solid line gives the best-fit power-law relation, with exponent ffflt = 1.933. 
The dashed line is the power law implied by the fit of our critical-O' model to 
the haloes' density and velocity-dispersion profiles in Figure|3] with e = 3 
and /?o = -0.10, ffciit = 1.96. Bottom panel: Relative residuals of (p/o'l) 
about the latter, critical-of power law assumed by our model fit. 
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Figure 5. Comparison of the best fit {blue soli d cun>e s) of equation 140bt 
to the high resolution ('billion particle') halo of lPiemand et ah 1 2005). For 
this fit, e = 3 was fixed a priori and only /?o, rg, and Mtot were allowed to 
vary. The vertical do t-dashed line in dicates the resolution linut of the sim- 
ulation, according to lPiemand et all The black dotted, green dash-dot, and 
red long-dashed curves correspond to the best fits of the profiles proposed 
by, respectively, Navarro et al. ( 1996) (see eq. 1431 ). Navarro et al- i 200i) 
(eq. (43), and lDiemand et alj 120031) (eq. l45l ). 



2, |log PraodC''/) - log pC^'i)!- In this case, the shape of the density 
profile at small r requires a slightly larger y6o than we found for Fig- 
ure|3| ySo - 0.03 (which, reasonably, is still nearly isotropic). This 
is one indication that dark-matter halo density profiles may not be 
exactly universal after all, even if the value of e is. Some of the 
scatt er found in y at the resolution limit of numerous simulations 
(e.g. iNavarro et aljEool iDiemand et"aLll2004bt iFukushige et alJ 
|2004) may in fact be real and, in our model at least, connected to 
non-universal halo kinematics. 

The top panel of Figure|5]plots the best fit of equation <40bt 
against the data for this one halo, now with the virial radius 
taken as the normalisation point; this model is given by the solid 
(blue) curve. For comparison with other fitting functions employed 
in the literature, we have also found the best-fitting NFW profile, 

P(r)/p2 = 4(r/r2)-\l + r/r2r\ (43) 

which we show as the dotted (black) curve. The best-fitting function 
of the type suggested by Navarro et al. ( 2004.) , 

p(r)/p2 = exp j- (2/a^) [(r/r2)''" - 1]) , (44) 

is shown as the dash-dot (green) curve. Here is a free parame- 
ter in the fit, which we find to be = .144. Finally, the fitting 
formula preferred bv lDiemand et alj i2005h . 

p(r)/p,, = 2'-^ (rlry \\ + r/r,)-' ** , (45) 

is drawn as the long-dashed (red) curve. Within the resolved ra- 
dial range of the simulation, r/rvh > 10"^^, there is little obvious 
difference between these fits in a plot of p vs. r. 

In the middle panel of Figure|5|we show the fractional residu- 
als from each of the four models fitted to p(r). The points joined by 



the solid blue curve denote the residuals from the fit of our model; 
residuals from the others just listed are in the same line types and 
colours as in the upper panel. The NFW model is clearly worse than 
any of the others, but the remaining three are competitive. The rms 
fractional density residual about the best NFW fit is 16%; that about 
the N avarro et aL function is 6.3%; that about the IPiemand et all 
(2005) formula is 5.2%; and that about our model with e = 3 and 
ySo = 0.03 is 6.0%. The main source of the slightly higher global 
scatter in our model vs. the fitting function of Diemand et al. is the 
single density point closest to the resolution limit of this simulation. 

It is noteworthy that our best-fit dynamical model is almost 
identical to the best-fitting Navarro et al. function within the re- 
solved radial range, as the latter has been shown to provide a very 
accurate description of many simulated d ark-matter haloes (e.g., 
iNavarro et all2004lDiemand et all2004d) . At some level, it is not 
surprising that either of these curves is an improvement over, say, 
the NFW profile, since the former both involve three free param- 
eters in p(r) (recall that we fixed e = 3 in fitting our model to 
this halo) while the NFW function contains only two. The advan- 
tage to our model, of course, is that the extra degree of freedom in 
fitting the density profile is also used simultaneously to predict the 
anisotropic kinematics of the halo. For example, from equation <31> 
we would predict a = o-cit = 1.938 in the p-cr relation js} for this 
halo. Perhaps coincidentally, this is almost exactl y the o-flt found 
in Figu relUfor the ten lower-resolution haloes from lPiemand et all 
teooia b). 

Finally, the bottom panel of Figure |5| shows the negative log- 
arithmic densit y slope estimated as a function of radius directly 
from the data of lPiemand et alJi2005h . against the behaviours pre- 
dicted by each of the four fits to the density profile. With 6 = 3 
and ySo = 0.03, our model has y — > 70 = 0.81 as r — > 0, and 



Dynamical insight into dark-matter haloes 1 1 



7 — > y„ = 3.44 as r — > oo. The rollover to 70 is very gradual how- 
ever (from eq. I36al . r] = 0.438 here) and at the resolution limit of 
the simulation we still find 7^1 for the fit. Apparently, still higher- 
resolution simulations are required to disting uish clearly bet ween 
our density model and others, such as that of lPiemand et all with 
dilferent asymptotic slopes in the limit r — » 0. 



5 SUMMARY AND DISCUSSION 

As a rule, investigations into the structure of simulated dark-matter 
haloes have focused on the nearly 'universal' shape of their den- 
sity profiles. This approach is plagued, however, by the fact that 
the haloes are only ever resolved over 2-3 decades in radius. Since 
the logarithmic density slope clearly varies smoothly over this en- 
tire range, in a way that is not anticipated from any first prin- 
ciples, it is then difficult to characterise unambiguously the true 
form of p(r) with necessarily ad hoc, empirical fitting functions. 
In particular, the value of the asymptotic central power-law slope 
(if a single, universal asymptote even exists in the limit r — » 0) 
remains ill-constrained. Thus, the starting point of our analysis 
in this paper was the remarkable, and much simpler, empirical 
fact iTavlor & N avarro 2001, see also Hansen 2004) that simulated 
haloes satisfy, over their whole resolved range, a p-cr relation of the 
general form 

P('-) -a 
oc r 

for cr, the radial component of velocity dispersion, and e - 3. (As 
it happens, using the total one-dimensional velocity dispersion also 
leads to a nearly power-law p-cr relation, with a slightly different 
value for a; see Taylor & Navarro, and Figure Q above.) A single 
power-law dependence such as this is clearly easier to recognise, 
and to quantify accurately, than is a radius-dependent curvature in 
the density profile. It is also the simplest nontrivial form that any 
halo relation could take, and hence the p-cr scaling may be the most 
fundamental aspect of dark-matter haloes. 

We should note here that, while p/cr' has the dimension 
of ph ase-space density, its interpretation as phase-space density 
{laylor & Navarro 2001; Williams et al. 2004), coarse-grained or 
not, is problematic*. On the other hand however, if the phase-space 
density were in some sense of power-law form, we might expect 
p/o"^ to closely, but not exactly, follow a power law, too. 

Our basic assumption has been that a power-law p-cr relation 
holds for all radii, not just those resolved by simulations. It is then 
possible to use this as a constraint on the spherical leans equation 
in order to derive the density profile of an equilibrium halo. In Sec- 
tion|5|we did just this, for a specialised case with e = 3 fixed and 
velocity isotropy assumed. This part of the analysis is thus similar 
to the origina l work of Taylor & Navarro 1 2001), and to aspects of 
IWiUiams etalj j2004 see also|Hanseii||2004j). However, our study 
has expanded considerably on these others, as we investigated the 
full space of all physical solutions to the isotropic Jeans equation 
under a power-law p-cr relation with arbitrary exponent a. Most 
importantly, we have found that only one solution p(r) exists which 
asymptotes to realistic power-law behaviour in the limits of both 
small and large radii. This solution occurs only for a single critical 

* Only if the functional form of the distribution of velocities were indepen- 
dent of radius, i.e. the distribution function sepai'able, would pju^ equal the 
average phase-space density at radius r. However, there is no a priori reason 
why this situation should be satisfied in dark-matter haloes. 



value of a and, for this one solution only, the density, velocity- 
dispersion, and enclosed-mass profiles are all fully analytic. 

In Section|3|we extended our analysis by allowing for arbitrary 
exponents e in the p-cr relation and, unlike any other detailled study 
before, for velocity anisotropy. We found that the problem remains 
analytically tractable for any e > 2, only if we adopt a p-cr relation 
involving the radial component cr,, of velocity dispersion rather 
than the total (empirically both are equally well motivated, see 
Fig.0 and if we assume the anisotropy parameter y6(r) = 1 -crg/cr; 
to be either spatially constant or linearly dependent on the logarith- 
mic density slope 7(r) = -din p/dln r. Under these reasonable 
assumptions, there always exists an exact analogue to the fully an- 
alytical density profile we found in the isotropic, e = 3 case. This 
solution again occurs for a particular, 'critical' a value. Moreover, 
for the realistic case of f} linearly increasing with 7 (corresponding 
to outwards ever more radially biased anisotropy), this is the only 
physical solution which resembles simulated dark-matter haloes. 

This result is gratifying for two reasons. First, it con- 
nects neatly with an en tirely independent, empirical finding by 
iHansen & Moord i2005h . who have argued on the basis of a vari- 
ety of simulations that indeed there exists a roughly linear relation 
between p{r) and 7(r) in collisionless haloes formed through nu- 
merous different processes. Second, it turns out that the slope of 
this linear 'fi-y relation' does not affect the critical value of a, nor 
does it affect the shape of the equilibrium p(r) profile itself. Rather, 
both of these things are influenced only by the central anisotropy 
ySo = P(r = 0). Specifically, for a p-cr relation with e = 3 (which 
ultimately does appear to be most appropriate), our solution of the 
Jeans equation with radially varying velocity anisotropy requires 



f- oc r-ii + ^A 
crj 

and has density 

7+10^ 



7+10^ / 2(2-ll„) yo 

p{r) oc X 9 1 -I- X ■> 
with x= r/ro, radial velocity dispersion 
cr^(r) OC II + I , 



and anisotropy 



P(r) = 1 



o-2(r) 



1 + ,^2(2-ft)/9 



(46a) 



(46b) 



(46c) 



(46d) 



for any y6o < 1, ySoo < 1. Thus, the halo density profile has a shallow 
power-law cusp, p — » ^-(7+io/?o)/9 jjj jjjg lirnjt r — > 0, and steep- 
ens monotonically but slowly to p — > ^-(3i-2ft)/'' 35 ^ _> 00. More 
general results for 6 5*= 3 in the p-cr relation, including analytical ex- 
pressions for the velocity-dispersion, enclosed-mass, and circular- 
velocity profiles, are given in equations <401 above (with auxiliary 
definitions in eqs. 1361 ). 

In Section |4| we fit our analytical, critical-cf density and 
velocity- dispersion profiles to thos e in ten dark-matter haloes sim- 
ulated bv iPiemand et al] i2004aibl) . When modelling these 'data', 
we assumed that all haloes obey a p-cr relation with a single value 
of e and have the same central velocity anisotropy /Jo, so that the 
halo density profile is necessarily universal. We were able to find a 
good fit to all ten haloes simultaneously for e = 3 and y6o = -0.10, 
only slightly different from isotropic; see Figure |3| For this com- 
bination of parameters, our model requires p/cr' oc r""* for self- 
consistency, and we showed in Figure |3] that this is a satisfactory 
description of the data. 
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In Section |4| we also fit our model to the density profile of 
an extremely h i gh-re solution ('billion-particle') halo simulated by 
IPiemand et alJ (2005), finding good agreement again with e = 3 
and a nearly isotropic ySo = 0.03. In this case the fitted density 
profile has a central power-law cusp p — > r^o.si ^s r — » 0, which 
steepens gradually to p — » as r — > oo; see Figure|5| At the res- 
olution limit of the Die mand et al.l simulation (about 0.001 of the 
virial radius), the fit has y - I, significantly larger than the asymp- 
totic cusp slope and still consistent with both the simulation data 
and the behaviour of other fitting functions. Within the resolved ra- 
dial range, our model fit is als o very c losely traced by the best fit of 
the ad hoc profile proposed bv lNavar ro et al. (2004). 

Our findings suggest a possible first-principles explanation for 
halo density profiles along the following line of arguments. The ini- 
tial distribution function before collapse was completely cold (a S- 
function in velocity space) and hence the phase-space density scale- 
free. The collapse and the subsequent process of violent relaxation 
(phase-space mixing) is driven by gravity alone, which cannot in- 
troduce any scale dependence. This implies that the phase-space 
density of the collapsed halo satisfies some form of scale invari- 
ance, suggesting that the ratio p(r)/<T^, which is closely related to 
the phase-space density, follows a power law (the general scale- 
invariant functional form). However, as our analysis has shown, if 
p(r)/(T^ of dark-matter haloes is any power of radius, then the con- 
dition of equilibrium (and physically sensible density profiles) re- 
quires the particular power law r"""" and, simultaneously, that the 
density profile follows that given above. 

From our fits to simulated haloes in Section |4] the /J-y rela- 
tion seems not as tight or universal as the p-tr relation. This is not 
very surprising, since different amounts of velocity anisotropy are 
required to stabilise different spatial halo shapes, resulting in some 
scatter between the /3-y relations of different haloes. Since only the 
central velocity anisotropy fig affects the value of a^^t, this scat- 
ter has little (or no) influence on the the p-cr relation as long as 
jSo is similar (or the same) for different haloes. This explains why 
haloes of different spatial shape still have very similar density pro- 
files. Whether /Jo is a universal parameter (close to zero) or whether 
there is some real scatter we cannot predict, but in Section|4|we ob- 
tained a good fit to ten different simulated haloes using a single 
value for Po (see Fig.|3j. 

Of course, our analysis is still restricted to spherical symmetry 
and ignores any halo substructure. However, substructure is unim- 
portant for the issue of the overall density profile, as simulations of 
dark-matter structure formation with suppressed small-scale power 
in the initial conditions still yield the same characteristic density 
profiles, but much less substructure (Moore et al. 1999). The issue 
of asphericity is more likely to be relevant and hence it is all the 
more remarkable that our spherical analysis gives such a good de- 
scription of the spherically averaged profiles. Evidently, anisotropy, 
which was ignored in previous studies, is presumably more impor- 
tant than asphericity, because the gravitational potential is always 
less aspherical than the density distribution. 
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